Wybór solwera w zależności od rodzaju problemu
Metoda elementów skończonych służyć może do rozwiązania wielu problemów inżynerskich. Ogólny schemat działania metody jest podobny, niezależnie od rozważanego problemu. W naszym przypadku dla uproszczenia wybieramy problem projekcji bitmapy.
Podstawowe etapy metody elementów skończonych wymienione są poniżej. Ostatnim etapem jest wybór algorytmu służącego do rozwiązania wygenerowanego układu równań. Algorytm ten zależy od struktury macierzy, która pośrednio zależy od rozważanego problemu inżynierskiego.
- Wybór równania różniczkowego opisującego problem inżynierski lub zjawisko fizyczne które chcemy symulować.
- Najprostszy możliwy problem inżynierski to problem projekcji, na przykład problem projekcji bitmapy, opisany w rozdziale pierwszym. W problemie tym równanie różniczkowe zastąpione jest problemem: Znaleźć funkcję \( u\in C^1(\Omega) \) taką która aproksymuje zadaną bitmapę, czyli \( u(x,y)\approx BITMAP(x,y),\textrm{ gdzie } \Omega \) to obszar na którym określony jest nasz problem inżynierski.
- Kolejny krok to przekształcenie naszego problemu do tak zwanej formy słabej, w której to nasze równanie \( u(x,y)\approx BITMAP(x,y),\textrm{ gdzie } \Omega \) próbkowane jest za pomocą wielu tak zwanych funkcji testujących \( v(x,y) \) i uśredniane za pomocą całek. W przykładowym problemie projekcji bitmapy, forma słaba ma następującą postać: \( \int_{\Omega} u(x,y)v(x,y)dxdy=\int_{\Omega} BITMAP(x,y)v(x,y)dxdy \) gdzie \( v(x,y) \) to funkcja testująca. Równanie to musi być spełnione dla wszystkich wybranych przez nas funkcji testujących.
- Następny etap to generacja siatki elementów skończonych pokrywających zadany obszar i wybranie funkcji bazowych rozpiętych na elementach skończonych. W etapie tym obszar na którym rozwiązujemy nasz problem inżynierski dzielony jest na elementy skończone. Następnie na elementach skończonych rozpinamy wybrane funkcje bazowe których kombinacja liniowa przybliżać będzie naszą poszukiwaną funkcję \( u \) W naszym przykładzie projekcji bitmapy opisanym w rozdziale pierwszym, obszar \( \Omega \) czyli prostokąt na którym rozpięta jest bitmapa, dzielony jest na \( N_x\times N_y \) elementów prostkątnych. Następnie na obszarze tym rozpinamy funkcje bazowe B-spline drugiego rzędu, używając knot vector'a [0 0 0 1 2 3 ... N_x-1 N_x N_x N_x] w kierunku osi \( x \) oraz knot vector'a [0 0 0 1 2 3 ... N_y-1 N_y N_y N_y] w kierunku osi \( y \). Mamy więc bazy dwuwymiarowych funkcji B-spline \( B_{i,j}(x,y)=B^x_i(x)B^y_j(y), i=1,...,N_x, j=1,...,N_y \).
- Dokonujemy teraz dyskretyzacji sformułowania słabego za pomocą wybranych funkcji bazowych rozpiętych na siatce obliczeniowej. Poszukiwaną funkcję \( u \) wyrażamy jako kombinacje liniową funkcji bazowych. Podobnie za poszczególne funkcje testujące wybieramy poszczególne funkcje bazowe. W przykładzie projekcji bitmapy opisanym w rozdziale pierwszym, funkcje bazowe B-spline wybrane w kroku poprzednim używane są teraz do aproksymacji rozwiązania \( u=\sum_{i=1,...,N_x;j=1,...,N_y} a_{i,j }B^x_i(x)B^y_j(y) \), podobnie jako funkcje testujące wybieramy funkcje B-spline \( B^x_k(x)B^y_l(y) \). Otrzymujemy w ten sposób następujący układ równań \( \begin{bmatrix} \int{B^x_{1,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int{B^x_{1,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int{B^x_{1,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \int{B^x_{2,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int{B^x_{2,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int{B^x_{2,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \vdots & \vdots & \vdots & \vdots \\ \int{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{1,p}B^y_{1,p}} & \int{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{N_x,p}B^y_{N_y,p}} \end{bmatrix}\begin{bmatrix} {\color{red}u_{1,1}} \\ {\color{red}u_{2,1}} \\ \vdots \\ {\color{red}u_{N_x,N_y}} \end{bmatrix} = \\ = \begin{bmatrix} \int BITMAP(x,y) {\color{blue}B^x_1(x)*B^y_1(y)} dx \\ \int BITMAP(x,y) {\color{blue}B^x_1(x)*B^y_2(y)} dx \\ \vdots \\ \int BITMAP(x,y) {\color{blue}B^x_{N_x}(x)*B^y_{N_y}(y)} dx \\ \end{bmatrix} \)
- Ostatni etap to rozwiązanie wygenerowanego układu równań.
Mamy sześć reprezentatywnych algorytmów solwerów służących do rozwiązywania układów równań wygenerowanych podczas obliczeń za pomocą metody elementów skończonych. Są to algorytm eliminacji Gaussa, algorytm LU faktoryzacji, algorytm solwera frontalnego, algorytm solwera wielo-frontalnego, algorytm solwera zmienno-kierunkowego, oraz algorytm iteracyjny.
Dla problemów dla których macierze są symetryczne i dodatnio określone, często używane są algorytmy iteracyjne bazujące na metodach Krylova, na przykład algorytm gradientów sprzężonych. Istnieją również algorytmy iteracyjne dla szerszej klasy problemów, np. algorytm GMRES.
Algorytmy solwerów iteracyjnych opisane zostały są w bardzo obszerny sposób w podręczniku Józefa Saada [1].
Liczba iteracji dla solwerów iteracyjnych zależy od rodzaju rozwiązywanego problemu, od struktury i własności macierzy po dyskretyzacji. W szczególności na zbieżność solwera typu gradientów sprzężonych wpływa współczynnik uwarunkowania macierzy, czyli iloczyn modułu z największej wartości własnej macierzy i modułu z największej wartości własnej macierzy odwrotnej. Jeśli współczynnik uwarunkowania macierzy jest za duży, wówczas tworzone są tak zwane preconditionery, czyli macierze przez które przemnaża się dodatkowo nasz układ równań w taki sposób żeby rozwiązanie było takie samo, ale współczynnik uwarunkowanie macierzy był lepszy. Niestety problem doboru preconditionera jest zaleźny od rozwiązywanego problemu (nie ma uniwersalnych zawsze skutecznych algorytmów).
Solwery dokładne z koleji potafią rozwiązać każdy problem, potrafią sfaktoryzować każdą macierz, jednak ich koszt jest zazwyczaj większy od kosztu solwerów iteracyjnych. Dla jednorodnych siatek dwuwymiarowych koszt ten jest rzędu \( {\cal O } ( N^{1.5 }) \). Dla jednorodnych siatek trójwymiarowych koszt ten jest rzędu \( {\cal O }(N^2) \). Wyjątkiem jest algorytm solwera zmienno-kierunkowego, którego koszt jest liniowy \( {\cal O}(N) \).
Algorytm solwera zmienno-kierunkowego może być zastosowany wtedy, kiedy macierz ma strukturę produktu Kroneckera. Dzieje się tak na przykład w problemie projekcji, lub podczas symulacji niestacjonarnych za pomocą metody explicite.
Algorytm solwera wielo-frontalnego, jest używany dla problemów dla których macierze nie są symetryczne lub nie są dodatnio określone. Algorytm solwera wielo-frontalnego jest szczególną implementacją algorytmu eliminacji Gaussa lub LU faktoryzacji.
Wersja iteracyjna algorytmu zmienno-kierunkowego, opisana w tym podręczniku używana jest do symulacji explicite na nieregularnych geometriach.